unit uMatrix;

interface

uses uNumber;

type

  TMatrixIndex = Byte;

  TMatrixBase = array of array of TNumber;

  TMatrix = class
  private
    FRowCount,
    FColumnCount: TMatrixIndex;
    FValueMatrix: TMatrixBase;
    FNumberClass: TNumberClass;
    FDeterminantSign: shortint;

    procedure CreateBase(var ABase: TMatrixBase; ARowCount, AColumnCount: TMatrixIndex);
    procedure CopyBase(ASource: TMatrixBase; var ATarget: TMatrixBase);
    procedure DestroyBase(AMatrix: TMatrixBase);

    function MinSize: TMatrixIndex;
    function GetStrValue(ARow, AColumn: TMatrixIndex): string;
    procedure SetStrValue(ARow, AColumn: TMatrixIndex; AValue: string);

    function GetRowString(ARowIndex: TMatrixIndex): string;

  public
    constructor Create(ARowCount, AColumnCount: TMatrixIndex; ANumberType: TNumberClass);
    destructor Destroy; override;

    procedure AddScalar(AScalar: TNumber); overload;
    procedure AddScalar(AScalarStr: string); overload;
    procedure SubtractScalar(AScalar: TNumber); overload;
    procedure SubtractScalar(AScalarStr: string); overload;
    procedure MultiplyByScalar(AScalar: TNumber); overload;
    procedure MultiplyByScalar(AScalarStr: string); overload;
    procedure DivideByScalar(AScalar: TNumber); overload;
    procedure DivideByScalar(AScalarStr: string); overload;

    procedure Add(AMatrix: TMatrix);
    procedure Subtract(AMatrix: TMatrix);
    procedure MultiplyBy(AMatrix: TMatrix);

    procedure AddToRow(ARow: TMatrixIndex; ANumber: TNumber); overload;
    procedure AddToRow(ARow: TMatrixIndex; ANumberStr: string); overload;
    procedure AddToColumn(AColumn: TMatrixIndex; ANumber: TNumber); overload;
    procedure AddToColumn(AColumn: TMatrixIndex; ANumberStr: string); overload;
    procedure MultiplyRow(ARow: TMatrixIndex; ANumber: TNumber); overload;
    procedure MultiplyRow(ARow: TMatrixIndex; ANumberStr: string); overload;
    procedure MultiplyColumn(AColumn: TMatrixIndex; ANumber: TNumber); overload;
    procedure MultiplyColumn(AColumn: TMatrixIndex; ANumberStr: string); overload;
    procedure DivideRow(ARow: TMatrixIndex; ANumber: TNumber); overload;
    procedure DivideRow(ARow: TMatrixIndex; ANumberStr: string); overload;
    procedure DivideColumn(AColumn: TMatrixIndex; ANumber: TNumber); overload;
    procedure DivideColumn(AColumn: TMatrixIndex; ANumberStr: string); overload;
    procedure AddRowToRow(ASourceRow, ATargetRow: TMatrixIndex);
    procedure AddColumnToColumn(ASourceColumn, ATargetColumn: TMatrixIndex);
    procedure AddMultipliedRowToRow(ASourceRow, ATargetRow: TMatrixIndex; AMultiplier: TNumber);
    procedure SwapRows(ARow1, ARow2: TMatrixIndex);
    procedure SwapColumns(AColumn1, AColumn2: TMatrixIndex);

    procedure Clear;
    procedure TopTriangle(ACanSwapColumns: Boolean = True);
    procedure BottomTriangle(ACanSwapColumns: Boolean = True);
    procedure Transpose;
    procedure Invert;
    procedure JordanEliminationMethod(ACanSwapColumns: Boolean = False);

    procedure CalculateDeterminant(ADeterminant: TNumber);

    procedure GetNumber(ARow, AColumn: TMatrixIndex; ANumber: TNumber);
    procedure SetNumber(ARow, AColumn: TMatrixIndex; ANumber: TNumber);

    property StrValue[ARow, AColumn: TMatrixIndex]: string read GetStrValue write SetStrValue;
    property RowStrings[ARowIndex: TMatrixIndex]: string read GetRowString;
    property RowCount: TMatrixIndex read FRowCount;
    property ColumnCount: TMatrixIndex read FColumnCount;
 end;

implementation

{ TMatrix }

procedure TMatrix.Add(AMatrix: TMatrix);
var
  r, c: TMatrixIndex;
begin
  //TODO: check size of AMatrix

  for r := 0 to FRowCount - 1 do
    for c := 0 to FColumnCount - 1 do
      FValueMatrix[r, c].Add(AMatrix.FValueMatrix[r, c]);
end;

procedure TMatrix.AddColumnToColumn(ASourceColumn, ATargetColumn: TMatrixIndex);
var
  r: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    FValueMatrix[r, ATargetColumn].Add(FValueMatrix[r, ASourceColumn]);
end;

procedure TMatrix.AddScalar(AScalarStr: string);
var
  scalar: TNumber;
begin
  scalar := FNumberClass.Create;
  try
    scalar.FromString(AScalarStr);
    AddScalar(scalar);
  finally
    scalar.Destroy;
  end;
end;

procedure TMatrix.AddMultipliedRowToRow(ASourceRow, ATargetRow: TMatrixIndex; AMultiplier: TNumber);
var
  c: TMatrixIndex;
  number: TNumber;
begin
  number := FNumberClass.Create;
  try
    for c := 0 to length(FValueMatrix[ATargetRow]) - 1 do
    begin
      number.CopyFrom(FValueMatrix[ASourceRow, c]);
      number.MultiplyBy(AMultiplier);
      FValueMatrix[ATargetRow, c].Add(number);
    end;
  finally
    number.Destroy;
  end;
end;

procedure TMatrix.AddRowToRow(ASourceRow, ATargetRow: TMatrixIndex);
var
  c: TMatrixIndex;
begin
  for c := 0 to FColumnCount - 1 do
    FValueMatrix[ATargetRow, c].Add(FValueMatrix[ASourceRow, c]);
end;

procedure TMatrix.AddScalar(AScalar: TNumber);
var
  r, c: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    for c := 0 to FColumnCount - 1 do
      FValueMatrix[r, c].Add(AScalar);
end;

procedure TMatrix.AddToColumn(AColumn: TMatrixIndex; ANumber: TNumber);
var
  r: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    FValueMatrix[r, AColumn].Add(ANumber);
end;

procedure TMatrix.AddToColumn(AColumn: TMatrixIndex; ANumberStr: string);
var
  number: TNumber;
begin
  number := FNumberClass.Create;
  try
    number.FromString(ANumberStr);
    AddToColumn(AColumn, number);
  finally
    number.Destroy;
  end;
end;

procedure TMatrix.AddToRow(ARow: TMatrixIndex; ANumberStr: string);
var
  number: TNumber;
begin
  number := FNumberClass.Create;
  try
    number.FromString(ANumberStr);
    AddToRow(ARow, number);
  finally
    number.Destroy;
  end;
end;

procedure TMatrix.AddToRow(ARow: TMatrixIndex; ANumber: TNumber);
var
  c: TMatrixIndex;
begin
  for c := 0 to FColumnCount - 1 do
    FValueMatrix[ARow, c].Add(ANumber);
end;

procedure TMatrix.BottomTriangle(ACanSwapColumns: Boolean);
begin
  //TODO: implement BottomTriangle
end;

procedure TMatrix.Clear;
var
  r, c: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    for c := 0 to FColumnCount - 1 do
      FValueMatrix[r, c].Clear;
end;

procedure TMatrix.CopyBase(ASource: TMatrixBase; var ATarget: TMatrixBase);
var
  r, c: TMatrixIndex;
begin
  for r := 0 to Length(ASource) - 1 do
  begin
    for c := 0 to Length(ASource[r]) - 1 do
    begin
      ATarget[r, c].CopyFrom(ASource[r, c]);
    end;
  end;
end;

constructor TMatrix.Create(ARowCount, AColumnCount: TMatrixIndex; ANumberType: TNumberClass);
begin
  inherited Create;
  FRowCount := ARowCount;
  FColumnCount := AColumnCount;
  FNumberClass := ANumberType;
  FDeterminantSign := 1;

  CreateBase(FValueMatrix, ARowCount, AColumnCount);
end;

procedure TMatrix.CreateBase(var ABase: TMatrixBase; ARowCount, AColumnCount: TMatrixIndex);
var
  r, c: TMatrixIndex;
begin
  SetLength(ABase, ARowCount, AColumnCount);
  for r := 0 to ARowCount - 1 do
    for c := 0 to AColumnCount - 1 do
      ABase[r, c] := FNumberClass.Create;
end;

destructor TMatrix.Destroy;
begin
  DestroyBase(FValueMatrix);
  inherited;
end;

procedure TMatrix.DestroyBase(AMatrix: TMatrixBase);
var
  r, c: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    for c := 0 to FColumnCount - 1 do
      AMatrix[r, c].Destroy;
  //AMatrix := nil;
end;

procedure TMatrix.DivideColumn(AColumn: TMatrixIndex; ANumber: TNumber);
var
  r: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    FValueMatrix[r, AColumn].DivideBy(ANumber);
end;

procedure TMatrix.DivideByScalar(AScalarStr: string);
var
  scalar: TNumber;
begin
  scalar := FNumberClass.Create;
  try
    scalar.FromString(AScalarStr);
    Self.DivideByScalar(scalar);
  finally
    scalar.Destroy;
  end;
end;

procedure TMatrix.DivideByScalar(AScalar: TNumber);
var
  r, c: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    for c := 0 to FColumnCount - 1 do
      FValueMatrix[r, c].DivideBy(AScalar);
end;

procedure TMatrix.DivideColumn(AColumn: TMatrixIndex; ANumberStr: string);
var
  number: TNumber;
begin
  number := FNumberClass.Create;
  try
    number.FromString(ANumberStr);
    DivideColumn(AColumn, number);
  finally
    number.Destroy;
  end;
end;

procedure TMatrix.DivideRow(ARow: TMatrixIndex; ANumber: TNumber);
var
  c: TMatrixIndex;
begin
  for c := 0 to FColumnCount - 1 do
    FValueMatrix[ARow, c].DivideBy(ANumber);
end;

procedure TMatrix.DivideRow(ARow: TMatrixIndex; ANumberStr: string);
var
  number: TNumber;
begin
  number := FNumberClass.Create;
  try
    number.FromString(ANumberStr);
    DivideRow(ARow, number);
  finally
    number.Destroy;
  end;
end;

procedure TMatrix.GetNumber(ARow, AColumn: TMatrixIndex; ANumber: TNumber);
begin
  ANumber.CopyFrom(FValueMatrix[ARow, AColumn]);
end;

function TMatrix.GetRowString(ARowIndex: TMatrixIndex): string;
var
  c: TMatrixIndex;
begin
  for c := 0 to FColumnCount - 1 do
  begin
    Result := Result + FValueMatrix[ARowIndex, c].ToString + ';';
  end;
end;

function TMatrix.GetStrValue(ARow, AColumn: TMatrixIndex): string;
begin
  Result := FValueMatrix[ARow, AColumn].ToString;
end;

procedure TMatrix.Invert;
var
  originalBase, extendedBase: TMatrixBase;
  r, c, i: TMatrixIndex;
begin
  if FRowCount <> ColumnCount then
  begin
    //TODO: throw an exception
  end;

  originalBase := FValueMatrix;

  // create extended base
  CreateBase(extendedBase, FRowCount, 2 * FColumnCount);
  CopyBase(originalBase, extendedBase);
  for i := 0 to FColumnCount - 1 do
  begin
    extendedBase[i, i + FColumnCount].BecomeMultiplicativeIdentity;
  end;

  FValueMatrix := extendedBase;
  FColumnCount := 2 * FColumnCount;

  try
    try
      JordanEliminationMethod(False);
    finally
      FValueMatrix := originalBase;
      FColumnCount := FColumnCount div 2;
    end;

    for r := 0 to FRowCount - 1 do
    begin
      for c := 0 to FColumnCount - 1 do
      begin
        FValueMatrix[r, c].CopyFrom(extendedBase[r, c + FColumnCount]);
      end;
    end;

  finally
    DestroyBase(extendedBase);
  end;
end;

procedure TMatrix.JordanEliminationMethod(ACanSwapColumns: Boolean);
var
  multiplier: TNumber;
  i, r, minSize: TMatrixIndex;
begin
  TopTriangle(ACanSwapColumns);

  //create neutral numbers (1) on the diagonal
  minSize := self.MinSize;
  multiplier := FNumberClass.Create;
  try

    for i := 0 to minSize - 1 do
    begin
      if not FValueMatrix[i, i].IsAdditiveIdentity then
      begin
        multiplier.CopyFrom(FValueMatrix[i, i]);
        multiplier.Invert;
        MultiplyRow(i, multiplier);
      end;
    end;

    //finish jordan method
    for i := minSize - 1 downto 1 do
    begin
      if not FValueMatrix[i, i].IsAdditiveIdentity then
      begin
        for r := i - 1 downto 0 do
        begin
          multiplier.CopyFrom(FValueMatrix[r, i]);
          multiplier.DivideBy(FValueMatrix[i, i]);
          multiplier.ChangeSign;
          AddMultipliedRowToRow(i, r, multiplier);
        end;
      end;
    end;

  finally

    multiplier.Destroy;

  end;
end;

function TMatrix.MinSize: TMatrixIndex;
begin
  if FRowCount > FColumnCount then
    Result := FColumnCount
  else
    Result := FRowCount;
end;

procedure TMatrix.MultiplyBy(AMatrix: TMatrix);
var
  newBaseRow, newBaseColumn, i: TMatrixIndex;
  newBase: TMatrixBase;
  product: TNumber;
begin
  if Self.RowCount <> AMatrix.ColumnCount then
  begin
    //TODO: throw an exception
  end;

  CreateBase(newBase, Self.RowCount, AMatrix.ColumnCount);
  product := FNumberClass.Create;

  try

    for newBaseRow := 0 to Self.RowCount - 1 do
    begin
      for newBaseColumn := 0 to AMatrix.ColumnCount - 1 do
      begin
        newBase[newBaseRow, newBaseColumn].Clear;
        for i := 0 to Self.ColumnCount - 1 do
        begin
          product.CopyFrom(FValueMatrix[newBaseRow, i]);
          product.MultiplyBy(AMatrix.FValueMatrix[i, newBaseColumn]);
          newBase[newBaseRow, newBaseColumn].Add(product);
        end;
      end;
    end;

    //replace old base with new base
    DestroyBase(FValueMatrix);
    FValueMatrix := newBase;
    FColumnCount := AMatrix.ColumnCount;

  finally
    newBase := nil;
    product.Free;
  end;
end;

procedure TMatrix.MultiplyByScalar(AScalar: TNumber);
var
  r, c: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    for c := 0 to FColumnCount - 1 do
      FValueMatrix[r, c].MultiplyBy(AScalar);
end;

procedure TMatrix.MultiplyByScalar(AScalarStr: string);
var
  scalar: TNumber;
begin
  scalar := FNumberClass.Create;
  try
    scalar.FromString(AScalarStr);
    MultiplyByScalar(scalar);
  finally
    scalar.Destroy;
  end;
end;

procedure TMatrix.MultiplyColumn(AColumn: TMatrixIndex; ANumberStr: string);
var
  number: TNumber;
begin
  number := FNumberClass.Create;
  try
    number.FromString(ANumberStr);
    MultiplyColumn(AColumn, number);
  finally
    number.Destroy;
  end;
end;

procedure TMatrix.MultiplyColumn(AColumn: TMatrixIndex; ANumber: TNumber);
var
  r: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    FValueMatrix[r, AColumn].MultiplyBy(ANumber);
end;

procedure TMatrix.MultiplyRow(ARow: TMatrixIndex; ANumberStr: string);
var
  number: TNumber;
begin
  number := FNumberClass.Create;
  try
    number.FromString(ANumberStr);
    MultiplyRow(ARow, number);
  finally
    number.Destroy;
  end;
end;

procedure TMatrix.MultiplyRow(ARow: TMatrixIndex; ANumber: TNumber);
var
  c: TMatrixIndex;
begin
  for c := 0 to FColumnCount - 1 do
    FValueMatrix[ARow, c].MultiplyBy(ANumber);
end;

procedure TMatrix.SetNumber(ARow, AColumn: TMatrixIndex; ANumber: TNumber);
begin
  FValueMatrix[ARow, AColumn].CopyFrom(ANumber);
end;

procedure TMatrix.SetStrValue(ARow, AColumn: TMatrixIndex; AValue: string);
begin
  FValueMatrix[ARow, AColumn].FromString(AValue);
end;

procedure TMatrix.CalculateDeterminant(ADeterminant: TNumber);
var
  i, minSize:  TMatrixIndex;
  originalMatrix: TMatrixBase;
begin
  //save the current matrix
  CreateBase(originalMatrix, FRowCount, FColumnCount);
  CopyBase(FValueMatrix, originalMatrix);

  FDeterminantSign := 1;
  TopTriangle;

  try
    ADeterminant.CopyFrom(FValueMatrix[0, 0]);

    //calculate determinant as multiplication of the diagonal numbers
    minSize := self.MinSize;
    for i := 1 to minSize - 1 do
      ADeterminant.MultiplyBy(FValueMatrix[i, i]);

    //if the determinant sign has changed
    if (FDeterminantSign = -1) then
      ADeterminant.ChangeSign;

  finally
    //restore original matrix (just pass the reference)
    DestroyBase(FValueMatrix);
    FValueMatrix := originalMatrix;
  end;
end;

procedure TMatrix.Subtract(AMatrix: TMatrix);
var
  r, c: TMatrixIndex;
begin
  //TODO: check size of AMatrix

  for r := 0 to FRowCount - 1 do
    for c := 0 to FColumnCount - 1 do
      FValueMatrix[r, c].Subtract(AMatrix.FValueMatrix[r, c]);
end;

procedure TMatrix.SubtractScalar(AScalar: TNumber);
var
  r, c: TMatrixIndex;
begin
  for r := 0 to FRowCount - 1 do
    for c := 0 to FColumnCount - 1 do
      FValueMatrix[r, c].Subtract(AScalar);
end;

procedure TMatrix.SubtractScalar(AScalarStr: string);
var
  scalar: TNumber;
begin
  scalar := FNumberClass.Create;
  try
    scalar.FromString(AScalarStr);
    SubtractScalar(scalar);
  finally
    scalar.Destroy;
  end;
end;

procedure TMatrix.SwapColumns(AColumn1, AColumn2: TMatrixIndex);
var
  r: TMatrixIndex;
  number: TNumber;
begin
  for r := 0 to FRowCount - 1 do
  begin
    number := FValueMatrix[r, AColumn1];
    FValueMatrix[r, AColumn1] := FValueMatrix[r, AColumn2];
    FValueMatrix[r, AColumn2] := number;
  end;
  FDeterminantSign := -1 * FDeterminantSign;
end;

procedure TMatrix.SwapRows(ARow1, ARow2: TMatrixIndex);
var
  c: TMatrixIndex;
  number: TNumber;
begin
  for c := 0 to FColumnCount - 1 do
  begin
    number := FValueMatrix[ARow1, c];
    FValueMatrix[ARow1, c] := FValueMatrix[ARow2, c];
    FValueMatrix[ARow2, c] := number;
  end;
  FDeterminantSign := -1 * FDeterminantSign;
end;

procedure TMatrix.TopTriangle(ACanSwapColumns: Boolean);
var
  r, c, i, minSize: TMatrixIndex;
  multiplier: TNumber;
  found: Boolean;
begin
  found := False;
  multiplier := FNumberClass.Create;

  try
  
    //iterate through diagonal
    minSize := self.MinSize;
    for i := 0 to minSize - 1 do
    begin
      //try to find neutral leading number (1) in lower rows
      if not FValueMatrix[i, i].IsMultiplicativeIdentity then
      begin
        r := i + 1;
        while r < FRowCount do
        begin
          if FValueMatrix[r, i].IsMultiplicativeIdentity then
          begin
            SwapRows(i, r);
            break;
          end;
          Inc(r);
        end;

        //if the neutral number wasn't found in the lower rows, try to find it in right cols
        if (not FValueMatrix[i, i].IsMultiplicativeIdentity) and (ACanSwapColumns) then
        begin
          c := i + 1;
          while c < FColumnCount do
          begin
            if FValueMatrix[i, c].IsMultiplicativeIdentity then
            begin
              SwapColumns(i, c);
              break;
            end;
            Inc(c);
          end;
        end;

      end;

      //try to find non zero leading number in lower rows (starting from the bottom)
      if FValueMatrix[i, i].IsAdditiveIdentity then
      begin
        found := False;
        r := FRowCount - 1;
        while r > i do
        begin
          if not FValueMatrix[r, i].IsAdditiveIdentity then
          begin
            SwapRows(i, r);
            break;
          end;
          Dec(r);
        end;
      end;

      //if non zero number still wasn't found, try to find it in other values
      if (FValueMatrix[i, i].IsAdditiveIdentity) and (ACanSwapColumns) then
      begin
        found := False;
        r := i;
        while (r < FRowCount) and (not found) do
        begin
          c := i + 1;
          while (c < FColumnCount) and (not found) do
          begin
            if not FValueMatrix[r, c].IsAdditiveIdentity then
            begin
              if (i <> r) then
                SwapRows(i, r);
              if (i <> c) then
                SwapColumns(i, c);
              found := True;
            end;
            Inc(c);
          end; //while c
          Inc(r);
        end; //while r
      end;

      //if all other values are zero
      if (FValueMatrix[i, i].IsAdditiveIdentity) then
        break;
 
      for r := i + 1 to FRowCount - 1 do
      begin
        //multiplier := aMatrix[r,i] / aMatrix[i,i]
        multiplier.CopyFrom(FValueMatrix[r, i]);
        multiplier.DivideBy(FValueMatrix[i, i]);

        //multiplier := -1*multiplier
        multiplier.ChangeSign;

        AddMultipliedRowToRow(i, r, multiplier);
      end;
    end;

  finally
    multiplier.Destroy;
  end;
end;

procedure TMatrix.Transpose;
var
  newMatrix: TMatrixBase;
  r, c: TMatrixIndex;
begin
  SetLength(newMatrix, FColumnCount, FRowCount);
  for r := 0 to FRowCount - 1 do
    for c := 0 to FColumnCount - 1 do
      newMatrix[c, r] := FValueMatrix[r, c];
  FValueMatrix := newMatrix;

  FRowCount := Length(newMatrix);
  FColumnCount := Length(newMatrix[0]);
end;

end.
